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Abstract 

The  project  was  initiated  to  try  and  obtain  answers  to  two  basic  questions: 

(1)  Do  non-equilibrium  energy  distributions  affect  macroscopic  reaction  rates? 

(2)  Can  these  effects  be  leveraged  for  control  in  design  of  scramjets? 

The  results  obtained  during  the  course  of  this  project  indicate  that  both  these  questions  can  be 
answered  affirmatively,  but  the  answers  were  surprising  and  non-intuitive.  Our  results  have 
shown  that  in  a  co-flow  hydrogen  jet  injector  vibrational  non-equilibrium  delays  ignition,  but  for 
a  highly  under-expanded  hydrogen  jet  in  a  supersonic  cross-flow  (as  in  the  Hyshot  configuration) 
vibrational  non-equilibrium  results  in  a  shorter  ignition  delay  with  the  flame  stabilizing  further 
upstream.  This  demonstrates  that  vibrational  non-equilibrium  can  be  used  to  improve  flame 
stability  in  a  suitably  designed  scramjet. 

A  high  speed  jet  flame  facility  was  designed  and  constructed.  Rayleigh  and  Raman  scattering 
measurements  were  made  in  this  flow  to  obtain  non-intrusive  temperature  measurements.  Using 
a  sheet  of  laser  light  single  shot  Rayleigh  scattering  was  used  to  obtain  instantaneous 
measurements  of  translational  temperature  in  a  plane.  The  Raman  scattering  measurements  were 
point  measurements  and  averaged  over  multiple  laser  shots  to  determine  independent  rotational 
and  vibrational  temperatures  using  high  dispersion  measurements  of  Stokes  Raman  scattering. 
The  Rayleigh  measurements  were  used  to  determine  the  probability  distribution  function  (pdf)  of 
the  turbulent  fluctuations.  Simulations  of  averaged  Raman  spectra  of  equilibrium  flows  using 
these  pdfs  were  used  to  quantify  the  spectral  distortion  induced  by  averaging.  The  distortions 
were  shown  to  be  small  provided  one  was  not  averaging  at  a  point  where  the  flame  is  only 
present  intermittently.  This  facility  was  used  to  show  that  in  the  mixing  layer  of  a  non-reacting 
supersonic  nitrogen  jet  in  air,  vibrational  non-equilibrium  of  nitrogen  in  the  mixing  layer  persists 
for  a  many  jet  diameters  downstream  of  injection.  Measurements  were  made  in  a  high-speed  jet 
flame  and  showed  clear  evidence  of  vibrational  non-equilibrium  upstream  of  the  flame,  and  even 
in  the  jet  core  downstream  of  the  flame  base.  Raman  measurements  on  both  nitrogen  and  oxygen 
gave  rotational  temperatures  that  were  in  excellent  agreement  with  each  other  and  with  the 
translational  temperature  measured  via  Rayleigh  scattering.  The  oxygen  vibrational  temperature 
was  found  to  be  in  equilibrium  with  the  rotational  and  translational  temperature,  but  the  nitrogen 
vibrational  temperature  was  out  of  equilibrium.  This  result  challenged  the  assumption  of  fast 
vibrational  energy  exchange  between  oxygen  and  nitrogen  that  was  commonly  used  in 
simulations. 
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A  highly  scalable  Quasi-Classical  Trajectory  (QCT)  code  was  developed  and  tested  using  a 
published  N4  potential.  Non-equilibrium  and  equilibrium  dissociation  rates  of  N2  were 
calculated.  The  equilibrium  rates  were  found  to  be  in  excellent  agreement  with  calculations  from 
another  group  using  the  same  potential  surface.  The  QCT  code  was  used  to  determine  non¬ 
equilibrium  reaction  rates  for  three  critical  reactions  in  the  hydrogen  combustion  mechanism. 
The  results  of  these  calculations  showed  that  the  chain-branching  combustion  reactions,  which 
are  exchange  reactions  rather  than  dissociation  reactions,  are  much  less  sensitive  to  vibrational 
temperature  than  is  assumed  in  the  standard  Park  two -temperature  model.  This  key  insight  is 
critical  to  understanding  why  vibrational  non-equilibrium  leads  to  enhanced  reaction  rates  in 
certain  flow  configurations.  The  relatively  elevated  translational  temperature  resulting  from  post¬ 
shock  gases  in  the  air-stream  being  vibrationally  cold  leads  to  a  net  increase  in  the  reaction  rates. 
The  hydrogen  fuel  in  the  under-expanded  injection  jet  is  vibrationally  “hot”,  but  this  does  not 
contribute  significantly  to  the  reaction  since  the  vibrational  temperatures  are  still  low.  The  QCT 
code  was  also  used  to  compute  V-V  exchange  probability  between  nitrogen  and  oxygen  at 
temperatures  relevant  to  our  laboratory  experiments.  The  results  confirmed  the  slow  V-V 
exchange  probability  that  produced  the  best  match  between  simulations  and  experiments.  These 
results  are  also  relevant  to  scramjet  inlet  and  isolator  conditions.  The  QCT  code  will  be  made 
available  to  the  community  via  a  web  repository. 

Detailed  numerical  simulations  were  performed  of  the  laboratory  jet  flame  configuration  and  of  a 
more  realistic  scramjet  combustor  configuration.  The  simulations  of  the  laboratory  flames 
showed  that  the  experimental  data  could  not  be  reproduced.  A  value  at  least  100  times  smaller 
than  the  standard  value  of  vibrational  exchange  probability  between  nitrogen  and  oxygen  was 
found  to  give  a  better  match  to  the  experimental  data.  This  value  was  then  confirmed  by  QCT 
calculations  of  nitrogen-oxygen  vibrational  energy  transfer  rates  over  a  wide  range  of 
temperatures.  Detailed  simulations  of  the  Hyshot  scramjet  configuration  using  representative 
flight  conditions  showed  that  vibrational  non-equilibrium  enabled  earlier  ignition  and  thus  a 
more  stable  flame  than  in  a  flow  with  complete  equilibrium. 
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Introduction 

This  report  summarizes  work  completed  under  this  grant  in  the  areas  of  experiment,  theory,  and 
numerical  flow  simulations.  In  this  final  report  we  focus  on  the  most  recent  results  obtained 
under  this  grant. 

1.  Experimental  effort 

We  designed  and  constructed  a  high  speed  jet  facility  to  enable  experimental  study  of  turbulent 
non-reacting  jets  and  turbulent  non-premixed  flames.  Laser  based  Rayleigh  and  high  dispersion 
Raman  scattering  were  used  to  probe  the  flow  non-intrusively. 

The  facility  has  a  small  axisymmetric  jet  in  a  high  temperature  air  coflow  and  is  mostly 
constructed  of  stainless  steel  to  prevent  rust  and  allow  for  high  temperature  operation.  The  jet 
and  coflow  issue  into  open  air,  and  the  coflow  velocity  is  less  than  1  m/s.  The  jet  gas  is  provided 
by  compressed  gas  cylinders  and  the  coflow  air  is  provided  by  a  blower  pulling  from  the  ambient 
room  air.  The  coflow  air  is  passed  through  a  filter  to  remove  the  dust  particles  larger  than  0. 1  pm 
in  order  to  facilitate  the  flow  diagnostics.  The  center  jet  flow  issues  through  a  CNC-machined 
stainless  nozzle,  which  can  be  replaced  to  provide  different  flow  velocities.  Figure  la  shows  the 
flow  facility  configuration  for  the  air  mixing  studies,  which  consists  of  a  Mach  1.5  round  jet, 
with  an  8  mm  exit  diameter,  and  centered  in  a  0.3  m  diameter  circular  coflow  chamber.  For  the 
reacting  flow  studies,  a  converging  nozzle  with  8  mm  exit  diameter  was  installed  and  a  conical 
shroud  reduced  the  coflow  diameter  to  100  mm.  The  reacting-flow  configuration  is  illustrated 
below  in  Fig.  lb.  This  configuration  enabled  the  collection  lens  used  in  Raman  and  Rayleigh 
scattering  measurements  to  be  brought  closer  to  the  probe  laser  beam  increasing  the  light 
collection  solid  angle  and  thus  improving  signal-to-noise  ratio  (SNR). 

This  facility  is  similar  to  those  used  in  previous  jet-in-coflow  studies1,2  except  that  the  coflow  is 
electrically  heated.  The  advantage  of  electrical  heating  is  that  it  does  not  alter  the  chemical 
composition  of  the  gases,  while  a  vitiated  coflow  would  include  combustion  products  that  could 
significantly  affect  the  vibrational  relaxation  process.  The  coflow  air  is  heated  by  a  pair  of  15 
kW  electric -resistance  flow-through  heaters.  The  maximum  temperature  of  the  heating  element  is 
rated  at  1200  K.  Two  type  K  thermocouples  inserted  radially  through  the  pipe  wall  75  mm  above 
the  heating  elements  measure  the  temperature  of  the  coflow  in  the  settling  chamber  and  allow 
feedback  control  of  the  heater  power.  The  jet  is  similarly  heated  by  a  pair  of  6  kW  electric 
heaters.  The  coflow  section  is  also  fitted  with  a  series  of  perforated  plates,  honeycomb  and 
screens  to  provide  flow  conditioning.  Translation  of  the  entire  assembly  is  provided  by  a  stepper 
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motor  and  screw-driven  translational  stage  for  simple  and  precise  (75  pm  accuracy)  scanning  of 
the  flow. 


Figure  1:  Cross-section  view  of  the  jet-in-coflow  facility,  (a)  Air- Air  mixing  configuration  with  Mach  1.5  nozzle 
installed,  and  (b)  reacting-flow  configuration  with  converging  nozzle  and  coflow  shroud. 

Air-Air  mixing  experiments 

In  one  set  of  experiments,  vibrational  non-equilibrium  was  studied  in  the  shear  layer  that 
developed  between  a  supersonic  air  jet  and  the  co-flow  of  heated  air.  The  operating  conditions 
are  summarized  in  Table  1.  The  air  jet  was  operated  in  a  perfectly  expanded  state  where  the  exit 
pressure  was  matched  to  the  ambient  pressure  of  the  room.  The  Raman  scattering  measurements 
were  made  in  the  axisymmetric  shear  layer  located  12  mm  downstream  of  the  jet  exit.  The  shear 
layer  thickness,  8,  shown  in  the  table  is  defined  by  the  5-95%  width  of  the  mean  temperature 
profile  (determined  by  Rayleigh  scattering).  The  characteristic  timescale  of  the  mixing  in  the 
shear  layer,  which  will  be  used  to  compare  with  the  relaxation  time  of  the  vibrational  energy, 
was  defined  as  the  eddy  turnover  time,  S/A U,  where  A  U  is  the  difference  between  the  speed  of 
the  jet  and  the  coflow.  This  timescale  is  very  short  relative  to  the  vibrational  relaxation  times  for 
N2  shown  below,  and  therefore  a  non-equilibrium  condition  is  expected. 


Table  1:  Operating  conditions  and  relevant  flow  parameters 
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Tjet 
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Air 
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210  K 
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0.51 

2.3  mm 

5  ps 

Spontaneous  Raman  scattering  setup 

Figure  2  shows  the  schematic  of  the  set  up  used  for  Raman  spectroscopy.  The  excitation  energy 
for  the  light  scattering  is  provided  by  a  frequency  doubled,  diode  pumped  Nd:YAG  laser 
(Coherent  Corona)  with  pulse  energy  of  4  mJ  at  532  nm  and  operated  at  10  kHz  repetition  rate. 
The  pulse  duration  is  -160  ns,  M2  is  25,  the  beam  width  is  5  mm  and  the  beam  divergence  is  5 
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mrad.  The  laser  light  is  switched  between  the  alignment  mode  and  the  experimental  modes  using 
the  combination  of  a  polarization  rotator  and  a  polarizing  beam  splitter.  In  the  experimental 
mode,  the  beam  splitter  is  removed  from  the  set  up  and  the  polarization  rotator  is  adjusted  to 
maximize  the  Raman  signal  in  the  direction  of  the  collection  optics.  The  laser  beam  is  focused 
using  a  40  cm  focal  length  lens.  The  scattered  light  is  collected  by  a  Nikon  105  mm  f/2.8  camera 
lens  and  focused  into  the  entrance  aperture  of  a  HoloSpec  f/1 .8  spectrograph  manufactured  by 
Kaiser  Optical  Systems,  Inc.  The  front  surface  of  the  collection  lens  is  positioned  about  22  cm 
from  the  nozzle  centerline.  The  camera  lens  used  is  corrected  for  chromatic  aberrations.  The 
elastically  scattered  light  is  filtered  out  by  a  SuperNotch  filter  with  O.D.  >  4.0  and  bandwidth  < 
350  cm-1.  Light  is  dispersed  by  the  HDG-608  transmission  grating,  which  has  an  average 
dispersion  of  -0.03  nm  per  pixel.  Both  the  notch  filter  and  the  grating  are  available  from  Kaiser 
Optical  Systems,  Inc.  The  signal  is  recorded  using  an  ICCD  camera  with  an  18  mm  wide  HBF 
Gen  III  intensifier.  The  intensifler’s  quantum  efficiency  is  44%  at  607  nm,  which  is  near  the  peak 
of  the  Q  branch  of  the  Stokes  signal  from  14N2.  The  CCD  array  has  1024  x  256  pixels  and  each 
pixel  has  dimensions  of  26  pm  x  26  pm.  The  spectral  coverage  is  595  to  620  nm,  the  slit  width 
of  the  spectrograph  is  100  pm  and  the  magnification  of  the  combined  collection  and  dispersion 
optics  is  approximately  0.8.  The  exposure  time  of  the  sensor  was  20  s,  i.e.,  each  data  set  is 
obtained  from  the  equivalent  of  800  J  of  incident  energy.  Such  high  incident  energy  is  necessary 
to  obtain  good  signal  to  noise  ratio  (SNR)  in  the  low-intensity  spectrally  resolved  O  and  S 
branches.  As  shown  in  Fig.  3,  the  intensity  of  the  O  and  the  S  branches  are  2-3  orders  of 
magnitude  less  than  that  of  the  Q  branch.  This  large  dynamic  range  is  made  possible  by  the  large 
vertical  sensor  area  on  which  the  signal  is  collected. 


Figure  2:  Schematic  of  the  Raman  set  up:  BD,  beam  dump;  BS,  beam  splitter;  PR,  polarization  rotator;  M,  mirror;  L, 

lens;  CL,  camera  lens;  and  NF,  notch  filter. 
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Wavelength  [nm] 

Figure  3:  Sample  Raman  spectrum  from  N2  in  the  jet  at  210  K.  Note  logarithmic  vertical  scale. 

The  spectral  calibration  procedure  is  as  follows.  The  raw  data  is  measured  in  terms  of  counts  as  a 
function  of  the  pixel  positions  along  the  dispersion  axis  of  the  CCD  array.  The  wavelength 
calibration  converts  pixels  to  wavelength,  and  the  intensity  calibration  converts  counts  to 
intensity.  The  wavelength  calibration  is  performed  by  comparing  the  pixel  positions  of  the 
emission  lines  from  a  neon  calibration  lamp  with  their  respective  positions  tabulated  in  the  NIST 
database.  The  pixel  positions  of  the  individual  neon  lines  are  obtained  to  sub  pixel  resolution  by 
curve  fitting  the  instrument  line  shape  model  to  the  recorded  spectrum.  The  instrument  line  shape 
function  is  modeled  by  convolving  a  trapezoid  function  and  a  Lorentzian.  The  free  parameters  of 
the  line  shape  function  are  the  half  base  width  and  the  half  top  width  of  the  trapezoid,  and  the 
Lorentz  width.  The  intensity  calibration  factors  are  obtained  by  taking  the  ratio  of  the  Planck’s 
function  at  1000  K  to  the  counts  recorded  from  a  blackbody  source  (Cl  Systems  SR-20)  at  the 
same  temperature.  Additional  details  on  the  data  calibration  procedure  are  described  in  Utsav,  et 
al.3 

Thermometric  Rayleigh  imaging  setup 

Planar  laser  Rayleigh  imaging  was  used  to  provide  an  instantaneous  measurement  of  the  gas 
temperature.  The  light  source  for  the  2D  Rayleigh  scattering  is  a  532  nm  frequency-doubled 
Nd:YAG  laser  (Spectra-Physics  GCR-150)  operated  at  10  Hz.  The  pulse  width  was  10  ns  and  the 
pulse  energy  was  300  mJ.  The  beam  was  first  expanded  using  a  pair  of  cylindrical  lenses  (with 
focal  lengths  -40  mm  and  130  mm,  respectively),  and  then  focused  with  a  500  mm  focal  length 
positive  cylindrical  lens.  The  resulting  laser  sheet  thickness  was  measured  by  traversing  a  knife- 
edge  through  the  beam  while  monitoring  the  power;  the  resulting  sheet  thickness  was  130  pm 
(FWHM).  The  sheet  height  was  10  ±  1  mm. 


Figure  4:  Schematic  diagram  of  the  Rayleigh  scattering  setup. 
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The  Rayleigh  scattering  was  imaged  using  a  back- illuminated  CCD  camera  (PixelVision 
SV512V1).  A  pair  of  Nikon  50mm  f/1 .4  lenses  were  used  in  a  “front-to-front”  configuration  and 
operated  at  infinite  conjugate  ratio  and  unity  magnification.  The  imaging  field  of  view  was  12 
mm  x  12  mm.  Resolution  of  the  imaging  system  was  characterized  by  measuring  the  step 
response  function  by  the  scanning  knife-edge  method.4  A  knife-edge  was  traversed  through  the 
object  plane  as  images  are  recorded  from  a  uniform  light  source.  By  plotting  the  response  in  a 
single  pixel  as  a  function  of  knife-edge  position,  the  line  spread  function  was  obtained.  Fitting  an 
error  function  to  the  response  curve  showed  the  full  width  at  half  maximum  to  be  29  pm.  The 
resulting  short  working  distance  of  the  lens  configuration  required  that  the  lenses  be  located 
within  the  hot  co flowing  air.  A  lens  protection  and  cooling  system  (not  pictured  in  Figure  4)  kept 
the  hot  co-flow  gas  from  damaging  the  lenses  and  camera.  Alteration  of  the  shear  layer  by  this 
intrusion  is  minimal  due  to  the  high  velocity  ratio  between  the  jet  and  coflow. 

A  set  of  sample  Rayleigh  scattering  images  is  shown  in  Fig.  5.  The  physical  region  being 
depicted  here  is  approximately  5  mm  x  6  mm  in  the  radial  and  streamwise  directions, 
respectively,  and  it  is  centered  12  mm  downstream  of  the  nozzle  exit.  The  denser  cold  jet  is  on 
the  left  while  the  hot  coflow  is  on  the  right.  The  shear  layer  exhibits  many  small  fingers  of  cold 
jet  fluid  reaching  into  the  hot  flow  without  the  coherent  “roller”  structures.  This  behavior  is 
consistent  with  previous  study  of  compressible  mixing  layer  structure.5,6 
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Figure  5:  Sample  Rayleigh  scattering  images 


Rayleigh  and  Raman  Temperature  Measurements 

Figure  6  shows  temperature  profiles  12  mm  downstream  of  the  nozzle  exit  obtained  using  the 
previously  described  techniques.  The  radial  coordinate  is  non-dimensionalized  in  a  manner 


consistent  with  previous  high-speed  mixing  studies:  rj  =  -  ^°5,  where  r0  5  is  the  radial  location 

of  the  mean  temperature  of  the  two  freestreams.  The  Rayleigh-based  profile  is  the  average 
temperature  field  computed  from  the  77  images  captured  during  the  3  minute  run.  The 
thermocouple  measurements  were  obtained  using  2  thermocouples  with  different  bead  diameters. 
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The  results  are  then  extrapolated  to  an  infinitesimally  small  bead  size  to  account  for  the  error  due 
to  bead  radiation.  Temperatures  from  the  Raman  scattering  measurements  are  shown  at  the 
discrete  locations  which  were  probed.  Vibrational  temperatures  are  omitted  where  the  air 
temperature  is  too  low  (<450  K)  for  a  meaningful  fit,  which  occurs  when  the  hot-band  intensity 
recedes  into  the  measurement  noise.  The  error  bars  in  Fig.  8  for  the  Raman  temperature 
measurements  were  determined  as  follows.  Residuals  between  the  collected  spectra  and  the  fit 
results  were  stored  in  a  library.  These  residuals  are  added  to  a  simulated  spectrum  at  known 
temperature  and  then  fit  for  temperature.  The  best  fit  value  of  temperature  from  the  spectral  fit  to 
the  noisy  simulated  spectrum  is  then  used  to  estimate  the  uncertainty  in  temperature  caused  by 
noise  in  the  data  and  temperature  variations  in  the  probe  volume  during  the  collection  time. 


Figure  6:  Various  temperature  measurements  as  a  function  of  radial  location  at  12  mm  downstream  of  the  nozzle 

exit. 

As  shown  in  Fig.  6,  there  is  good  agreement  between  the  Raman,  Rayleigh  and  thermocouple 
temperature  measurements  outside  the  shear  layer.  There  is  also  minimal  difference  between  7> 
and  Tr  away  from  the  shear  layer,  which  is  expected  for  flow  that  is  in  equilibrium.  On  the  hot 
side  of  the  shear  layer,  the  effect  of  slow  vibrational  relaxation  is  clearly  evident.  The  vibrational 
temperature  within  the  shear  layer  remains  higher  than  the  flow  temperature  as  the  energy 
relaxes  more  slowly  from  the  initial  inflow  step  function  profile.  This  can  be  imagined  as  the 
vibrational  energy  “lagging”  behind  the  rapid  cooling  of  the  flow  by  remaining  nearer  to  the 
freestream  temperature.  The  rotational  temperature  inferred  from  the  Raman  scattering 
measurements  agrees  with  the  Rayleigh-derived  temperature  to  within  50  K  at  all  locations.  One 
should  note  that  the  Raman  rotational  temperature  measurements  are  consistently  lower  than  the 
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Rayleigh  scattering  results  in  the  shear  layer.  This  may  not  be  an  error,  but  instead  could  be  a 
result  of  the  space  and  time  averaging  process.  Because  higher  density  fluid  scatters  light  more 
effectively,  the  collected  Raman  spectra  will  be  weighted  toward  the  colder  side  of  the 
temperature  probability  density  function  (PDF)  within  the  measurement  volume.  The  effect  of 
this  spatiotemporal  averaging  is  evaluated  further  in  the  following  section. 

Distinction  between  Compound  and  Non- equilibrium  Spectra 

One  drawback  of  using  time-averaged  Raman  spectroscopy  in  a  turbulent  flow  is  the  possibility 
of  inferring  false  levels  of  non-equilibrium  owing  to  the  effect  of  averaging  on  the  spectra.  In 
particular,  since  the  Raman  spectral  intensities  are  a  non-linear  function  of  temperature,  then  the 
average  spectrum  in  a  flow  with  fluctuating  temperature,  is  not  the  same  as  the  spectrum 
associated  with  the  average  temperature.  It  is  critical  that  the  magnitude  of  this  effect  is 
determined  before  conclusions  can  be  drawn  regarding  the  degree  of  non-equilibrium.  Our 
approach  is  to  create  a  compound  spectrum,  which  is  the  average  of  equilibrium  spectra  at 
different  temperatures,  and  then  to  fit  the  compound  spectrum  to  obtain  estimates  of  rotational 
and  vibrational  temperature.  If  the  non-linear  weighting  effect  is  negligible  then  the  fitting 
process  will  show  no  non-equilibrium  in  the  compound  spectrum. 

The  compound  equilibrium  spectrum  is  simulated  by  adding  the  spectra  computed  at 
temperatures  corresponding  to  the  actual  temperature  fluctuations  in  the  shear  layer  as  measured 
by  Rayleigh  scattering.  The  Rayleigh  temperature  measurements,  with  higher  spatial  and 
temporal  resolution,  provide  information  on  both  the  spatial  and  temporal  temperature 
distributions  within  a  given  Raman  probe  volume.  As  stated  previously,  the  Rayleigh  scattering 
measurements  have  a  spatial  resolution  of  29  pm  in  the  imaging  plane  and  an  exposure  equal  to 
one  laser  pulse  width,  10  ns.  Conversely,  the  Raman  spectra  are  recorded  from  a  1.1  x  0.33  mm 
region  and  integrated  over  20  s.  The  probability  density  functions  of  temperature  for  several  of 
the  Raman  sampling  locations  are  shown  in  Fig.  7.  These  sampling  regions  simulate  the  physical 
region  from  which  spectra  are  collected  by  the  spectrograph.  Equilibrium  Raman  spectra  are 
calculated  at  each  of  the  temperatures  shown  in  the  PDF  in  Fig.  7  and  are  averaged  to  generate  a 
compound  spectrum  at  each  location.  In  order  to  make  the  simulated  spectrum  more  similar  to 
the  experimental  data,  the  residuals  from  a  curve  fit  to  the  Raman  data  are  added  to  the  simulated 
compound  spectrum.  The  non-equilibrium  Raman  code  is  then  used  to  fit  the  simulated  noisy 
compound  spectrum.  We  find  that  the  non-linear  weighting  does  indeed  cause  some  apparent 
non-equilibrium  since  Tr  and  T y  differ;  however,  this  effect  is  less  than  50  K  at  all  locations. 
This  difference  is  significantly  less  than  the  difference  between  these  values  obtained  from  the 
experimental  data  shown  in  Fig.  6,  which  indicates  clearly  that  the  recorded  spectrum  is  not  just 
the  superposition  of  equilibrium  Raman  spectra. 
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Figure  7:  Rayleigh  temperature  PDFs  at  various  Raman  probe  regions  with  the  temperatures  obtained  from  spectral 

fitting  the  average  Raman  spectra  over  the  PDF. 

It  is  worth  noting  that  the  compressible  nature  of  the  mixing  layer  and  subsequent  suppression  of 
large-scale  “roller”  structures  may  help  to  reduce  the  temperature  fluctuations  in  the  shear  layer 
and  hence  the  magnitude  of  the  non-linear  weighting  error.  The  narrowing  of  local  PDFs  with 
increasing  compressibility  has  been  found  in  both  planar7,8,9  and  axisymmetric  mixing  layers.10 
In  other  flows  where  the  fluctuations  may  be  higher,  the  local  PDF  for  each  measurement  may  be 
broad  enough  that  the  difference  between  7>  and  Tr  in  a  compound  spectrum  may  be  of  the  same 
order  as  the  fitting  results  from  the  actual  Raman  scattering  data.  Because  of  this,  the  local 
temperature  fluctuations  must  be  quantified  when  extending  this  technique  to  other  flow 
environments. 

Effect  of  CO2  Addition 

Figure  8  shows  the  result  of  several  spectra  rapidly  collected  at  one  location  near  the  center  of 
the  shear  layer.  Samples  were  taken  first  with  the  air-air  mixing  case  presented  above  and  then 
CO2  was  seeded  into  both  streams  at  a  mole  fraction  of  approximately  5  and  10%.  The 
suppression  of  the  hot  band  at  606.5  nm,  as  a  result  of  the  addition  of  CO2,  signifies  the 
reduction  in  the  v  =  1  population.  Spectral  fits  showed  that  the  difference  between  7>and  Tr 
decreased  from  135  K,  without  any  addition  of  CO2,  to  90  K  with  the  addition  of  10  %  CO2.  This 
reduction  in  non-equilibrium  confirms  that  the  presence  of  CO2  accelerated  the  V-T  relaxation  of 
N2  as  expected,  while  also  reinforcing  the  fact  that  true  V-T  non-equilibrium  is  being  measured  in 
this  flow. 
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Figure  8:  Effect  of  C02  addition  on  vibrational  hot  band  population 
Competing  Excitation  and  De-excitation  Processes 

Thermal  energy  transfer  between  the  two  streams  occurs  as  hot  and  cold  packets  of  fluid  in  local 
vibrational-translational  equilibrium  are  entrained  from  the  low  and  high  speed  sides  of  the  layer, 
respectively.  If  one  assumes  this  process  to  have  equal  forward  and  reverse  relaxation  rates,  the 
temperature  profiles  would  be  expected  to  evolve  in  a  manner  similar  to  a  one-dimensional 
diffusion  process.  In  this  scenario,  since  density  effects  have  been  ignored,  the  process  appears  to 
behave  in  a  symmetrical  fashion.  This  would  be  analogous  to  a  shear  layer  in  which  the  mass 
entrainment  from  both  sides  of  the  layer  is  equal  so  that  the  mean  temperature  inside  the  layer 
would  match  the  mean  of  the  two  freestreams.  The  entrainment  ratio  of  the  shear  layer  can  be 
predicted  by  the  relation  given  by  Dimotakis.11  Because  the  velocity  ratio  (r)  in  our  flow  is  very 
low,  the  relation  simplifies  to 

Era(r^0)  =  1.68/.s',/2, 

where  5  is  the  density  ratio  across  the  shear  layer.  We  estimate  the  mass  entrainment  ratio,  Em, 
for  the  air-air  mixing  case  to  be  approximately  3.3,  which  indicates  a  strong  bias  in  mass  flux 
entering  the  shear  layer  from  the  high-speed  (cold)  side.  This  entrainment  asymmetry  will  reduce 
the  mean  temperature  of  the  shear  layer  so  that  the  hot  side  will  now  be  much  further  from  the 
equilibrium  temperature.  This  effect  is  illustrated  in  the  cartoon  profiles  in  Fig.  9,  which 
illustrate  how  the  apparent  non-equilibrium  would  appear  larger  in  magnitude  for  the  cooling 
process  for  entrainment  ratios  greater  than  one,  which  agrees  with  the  temperature  profiles 
presented  above. 
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Cross-stream  coordinate 


Figure  9:  Theoretical  depiction  of  effect  non-unity  entrainment  ratio  has  on  temperature  profiles. 


Within  each  probe  volume,  there  is  a  local  averaging  effect  as  fluid  packets  that  originate  on 
either  side  of  the  layer  are  present.  This  averaging  effect  can  reduce  the  apparent  magnitude  of 
non-equilibrium  as  molecules  with  excess  vibrational  energy  and  those  with  a  vibrational  deficit 
will  both  be  present  and  mixed  at  a  molecular  level.  To  isolate  this  effect,  dissimilar  species  can 
be  used  in  order  to  track  molecules  from  each  side  of  the  layer.  Because  the  jet  temperature  is  too 
low  to  probe  with  the  current  technique,  it  is  convenient  to  continue  to  use  air  on  the  hot  side  of 
the  layer  and  introduce  the  Raman  inactive  species,  argon,  in  the  jet.  This  allows  the  tracking  of 
vibrationally  hot  N2  molecules  as  they  are  entrained  and  mixed  with  the  cold  jet. 


Figure  10:  Temperature  profiles  for  the  air-air  mixing  and  argon-air  mixing  cases. 

In  Fig.  10,  the  results  are  presented  for  an  argon-air  case  and  compared  with  the  air-air  case 
shown  in  Fig.  6.  The  freestream  conditions  are  also  presented  in  Table  2.  The  conditions  for  the 
argon-jet  case  are  different  from  that  of  the  air -jet  case  since  the  ratio  of  specific  heats  is 
different  and  so  the  nozzle  must  be  run  off-design  to  achieve  pressure  matched  conditions  at  the 
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exit.  The  profile  is  truncated  where  the  local  concentration  of  nitrogen  is  too  low  for  the  fitting 
code  to  converge  to  a  meaningful  result. 

Table  2:  Experimental  conditions  during  the  air-air  and  argon-air  mixing  cases 


Air-Air 

Argon- Air 

Tjet 

210  K 

170  K 

T co flow 

850  K 

850  K 

Mjet 

1.50 

1.55 

Jjet 

1.4 

1.67 

Red 

8xl05 

6.3xl05 

The  non-equilibrium  can  be  observed  deeper  into  the  layer  when  vibrationally  active  molecules 
only  originate  on  the  hot  side  of  the  layer.  This  confirms  that  the  hot-to-cold  relaxation  process 
actually  produces  a  stronger  non-equilibrium  than  is  suggested  by  Fig.  6.  Therefore,  a  converse 
process  of  vibrationally  cold  molecules  relaxing  to  a  higher  temperature  must  also  exist  in  the 
layer.  Because  this  is  an  effect  that  is  present  at  the  molecular  level,  and  the  molecules  are 
distributed  homogeneously  within  the  probe  volume,  one  cannot  resolve  this  by  improving  the 
spatial  resolution  of  the  measurement.  In  a  shear  layer  with  identical  species  on  each  side,  the 
only  true  way  to  separate  the  two  processes  would  be  to  track  on  which  side  of  the  layer  each 
molecule  originated  from.  Because  this  is  not  possible  in  a  real  flow,  this  averaging  effect  will  be 
present  in  any  measurement  made  in  a  layer  with  vibrationally  active  species  in  both  freestreams. 

Summary  of  jet  mixing  studies 

The  effect  of  slow  vibrational  relaxation  in  a  high-speed  shear  layer  was  measured  using  a 
combination  of  time  averaged  spontaneous  Raman  scattering  and  instantaneous  planar  Rayleigh 
scattering.  The  spectral  fitting  Raman  procedure  is  sufficient  to  detect  non-equilibrium  when  the 
local  temperature  is  above  450  K.  It  was  only  possible  to  measure  the  non-equilibrium  between 
the  vibrational  and  the  rotational  energy.  It  was  not  possible  to  identify  non-equilibrium 
vibrational  distributions  because  only  the  fundamental  and  the  first  hot  bands  were  significantly 
populated.  The  introduction  of  CO2,  a  quenching  species,  suppressed  the  hot  band  as  expected, 
and  thus  provided  validation  that  the  non-equilibrium  observed  in  the  air-air  case  is  physical 
rather  than  an  artifact  of  the  measurement  technique.  The  overall  shape  of  the  vibrational  and 
rotational  temperature  profiles  was  shown  to  be  a  result  of  two  competing  non-equilibrium 
processes  occurring  within  the  shear  layer.  This  averaging  effect,  which  occurs  on  a  molecular 
level,  will  be  present  in  any  shear- layer- induced  non-equilibrium  process  that  has  the  same 
vibrationally  participating  species  on  both  sides  of  the  layer. 

High-speed  jet  flame  experiments 

In  the  following  study,  the  same  techniques  were  applied  to  high-speed  jets  of  combustible  gas 
mixing  with  hot  air  coflow.  Several  flame  conditions  were  run,  which  are  summarized  in  Table  3 
below.  For  each  run,  the  coflow  temperature  was  kept  at  1000  K.  The  jet  velocity  values  in  Table 
3are  calculated  based  on  input  mass  flow  rates  and  thermocouple  measurements.  The  flow 
conditions  in  this  study  were  chosen  such  that  the  flame  stabilized  at  approximately  the  same 
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liftoff  height  from  the  nozzle  exit.  Each  condition  provides  a  flame  that  auto-ignited  immediately 
at  the  onset  of  fuel  flow.  The  first  two  fuel  mixtures  listed  consist  of  H2  with  diluents  N2  and  Ar 
to  simplify  the  non-equilibrium  relaxation  mechanism  to  a  system  of  diatomic  and  atomic 
species,  minimizing  the  computational  expense  of  modeling  these  cases  in  companion  CFD 
simulations.  The  H2/N2  and  H2/Ar  cases  also  provide  a  comparison  pair  with  and  without  the 
presence  of  vibrationally  cold  N2  from  the  jet.  The  H2/CH4  fuel  mixture  is  the  same  as  in  [69] 
and  provides  a  Rayleigh  cross  section  which  is  identical  to  that  of  air.  This  cross  section 
matching  allows  for  the  application  of  Rayleigh  thermometry  to  quantify  the  effect  of  temporal 
averaging  in  the  Raman  measurements  in  a  similar  fashion  to  the  air-air  mixing  study.  Figure  1 1 
below  shows  the  average  lifted  flame  location  for  the  H2/N2  flame  case  as  marked  by  a  long 
exposure  OH*  chemiluminescence  image.  The  flame  lifts  approximately  5  diameters  from  the 
nozzle  exit  in  this  case  and  is  highly  turbulent.  Owing  to  the  much  higher  speed  of  sound  for  the 
combustible  mixtures,  all  the  cases  for  this  study  were  subsonic.  For  all  cases  the  stoichiometric 
mixture  fraction,  Zstoich,  is  relatively  small  and  so  the  flame  resides  mainly  on  the  low-speed  side 
of  the  shear  layer. 
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Figure  11:  Time-average  image  (false  color)  of  OH*  chemiluminescence  for  the  H2/N2  flame  case 


Table  3:  Summary  of  run  conditions 


Case 

Jet  (by  mol.) 

Co  flow 

Ujet  [m/s] 

Red 

Tojed  K] 

Zstoich 

H2/N2 

68%  H2,  32%  N2 

Air 

540 

75,100 

560 

0.181 

H2/Ar 

75%  H2,  25%  Ar 

Air 

505 

77,000 

570 

0.183 

H2/CH4 

62%  H2,  38%  CH4 

Air 

300 

27,400 

550 

0.048 
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For  this  study,  several  changes  were  made  to  the  Raman  and  Rayleigh  scattering  measurements. 
The  laser  used  was  a  frequency  doubled  diode-pumped  Nd:YLF  laser  (Coherent  Evolution  90) 
operated  at  5  kHz  repetition  rate  and  an  average  power  of  60W.  The  105mm  f/2.8  camera  lens 
was  no  longer  used  for  the  collection  optics.  Instead,  a  pair  of  172.2  150  mm  diameter  achromatic 
lenses  (Special  Optics  54-120-260)  operated  at  infinite  conjugate  ratio  was  used.  The  smaller  f- 
number  enabled  reduced  integration  times  to  achieve  the  same  signal  levels  as  with  the  previous 
collection  optics.  For  this  study,  we  also  obtained  a  grating  which  allowed  for  the  measurement 
of  O2  Stokes  Raman  scattering  by  giving  spectral  coverage  of  562  to  583  nm.  Owing  to  the  heat 
generated  by  the  flame,  the  paired  lens  setup  was  not  used  for  the  Rayleigh  scattering 
measurements  obtained  in  this  study.  Instead  a  forward-facing  Nikon  50mm  171.2  lens  was  used 
with  20  mm  of  extension  to  obtain  sufficient  magnification  while  maintaining  a  large  enough 
working  distance  to  protect  the  lens  and  camera  from  thermal  damage.  The  incident  light  for  the 
Rayleigh  scattering  was  provided  by  frequency-doubled  Nd:YAG  laser  delivering  1 J  per  pulse  at 
a  10  Hz  repetition  rate.  The  scattered  light  was  imaged  by  a  CCD  camera  (PCO  1400)  with  2x2 
pixel  binning  to  increase  the  signal-to-noise  ratio.  The  field  of  view  for  this  setup  was 
approximately  17  mm  and  10  mm  in  the  radial  and  axial  directions  respectively. 

Confirmation  of  non-equilibrium  in  measured  spectra 

In  the  air-air  mixing  study,  the  Raman-based  measurements  of  non-equilibrium  were  aided  by 
the  use  of  Rayleigh  thermometry  that  was  enabled  by  the  constant  Rayleigh  cross-section  of  the 
air-air  system.  The  same  analysis  can  be  applied  to  the  current  H2/CH4  Rayleigh  thermometry 
results  to  confirm  that  temperature  fluctuations  within  the  probe  volume  do  not  lead  to  excessive 
false  non-equilibrium.  Temperature  values  for  regions  of  the  flow  field  of  the  same  physical 
dimension  as  the  Raman  measurement  were  extracted  from  a  set  of  1000  Rayleigh  images.  The 
probability  density  function  within  each  of  these  regions  is  shown  in  Fig.  12.  Equilibrium  Raman 
spectra  are  calculated  and  summed  to  simulate  the  spatial  and  temporal  averaging  that  occurs 
during  a  Raman  measurement.  The  resultant  composite  spectra  are  then  fitted  for  vibrational  and 
rotational  temperatures,  which  are  also  presented  for  each  location  in  Fig.  12.  This  modeling 
exercise  shows  that  the  actual  temperature  fluctuations  cause  a  maximum  discrepancy  between 
Tv  and  Tr  of  at  most  15  K,  which  is  significantly  less  than  the  magnitude  of  non-equilibrium  seen 
in  the  collected  spectra;  however,  this  is  not  expected  to  be  the  case  in  a  region  over  which  the 
flame  is  only  present  intermittently.  Figure  13  shows  sample  Rayleigh  thermometry  images  from 
the  H2/CH4  flame,  which  illustrate  the  highly  intermittent  nature  of  the  flame.  The  flame  liftoff 
height  fluctuates  over  a  distance  that  is  larger  than  the  imaging  field  of  view.  A  time-averaged 
Raman  measurement  made  in  this  region  would  be  integrating  over  gas  temperature  variations  of 
over  1000  K  and  could  not  be  used  reliably  to  detect  vibrational  non-equilibrium. 
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Figure  12:  Temperature  probability  density  functions  extracted  from  Rayleigh  thermometry  images  along  with 
temperatures  obtained  by  applying  the  fitting  code  to  similarly  distributed  equilibrium  spectra.  The  profde  shown  is 

6  diameters  downstream  of  nozzle  exit. 


Raman  measurements  in  hydrogen  flames 

Using  hydrogen  as  fuel  simplifies  the  combustion  chemistry  and  allows  for  simpler  models  to  be 
used  for  CFD  validation.  The  lack  of  carbon-containing  polyatomic  species  also  significantly 
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simplifies  the  modeling  of  V-V  and  V-T  energy  transfer.  Radial  shear  layer  profiles  were  taken 
at  various  downstream  locations  (Fig.  14)  to  characterize  the  non-equilibrium  in  the  shear  layer 
upstream  of  the  flame. 


—  H2  measurements 

—  CH4  measurements 


Figure  14:  Schematic  diagram  of  the  Raman  scattering  measurement  locations  (to  scale). 


Figure  15a  shows  vibrational  and  rotational  temperature  measurements  at  a  location  that  is  in  the 
non-reacting  mixing  region  upstream  of  the  lowest  extent  of  the  unsteady  flame  base.  The 
rotational  temperatures  of  N2  and  O2  agree  very  well,  because  R-T  relaxation  is  very  fast.  This 
consistency  between  measurements,  which  were  performed  with  separate  runs  and  using  a 
different  spectral  region,  provides  verification  of  the  measurement  accuracy.  From  the 
vibrational  temperature  measurements  in  Fig.  15a  it  is  apparent  that  the  N2  is  out  of  thermal 
equilibrium  within  the  shear  layer,  as  seen  in  our  previous  measurements,12  but  interestingly  the 
O2  seems  to  be  in  complete  equilibrium  throughout  the  mixing  region.  This  difference  among  N2 
and  O2  would  suggest  that  the  V-V  intermolecular  transfer  between  these  two  species  is  very 
inefficient,  which  would  lead  to  each  species  relaxing  at  its  own  V-T  rate.  This  conclusion  agrees 
with  the  study  of  Cutler  et  al.,13  who  found  that  Tv^2  and  7fo2  were  quite  different  in  the  non¬ 
equilibrium  flow  exiting  a  Mach  2  heated  wind  tunnel  nozzle. 

The  error  bars  in  the  figure  indicate  the  precision  uncertainty  (95%  confidence)  that  was 
determined  using  repeated  measurements.  Measurements  are  presented  for  radial  locations  where 
the  signal-to-noise  ratios  are  high  enough  that  accurate  fits  could  be  made.  In  the  near  field, 
where  there  is  still  a  potential  core,  there  is  no  O2  in  the  jet  fluid,  and  so  the  O2  Raman  signals 
diminish  faster  than  the  N2  signals  as  the  probe  volume  is  translated  toward  the  jet  axis. 
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Figure  15:  a)  Radial  trace  of  N2  and  02  temperatures  with  H2/N2  fuel  and  b)  radial  trace  of  N2  temperatures  with 
H2/N2  and  H2/Ar  fuel.  Measurements  were  all  taken  at  an  axial  distance  of  2  jet  diameters  downstream  of  the  nozzle 

in  the  pre-flame  region. 

As  with  the  air-air  mixing  study,  the  apparent  non-equilibrium  strength  in  the  H2/N2  case  is 
reduced  by  mixing.  As  previously,  the  N2  in  the  jet  is  replaced  with  argon  to  isolate  the  hot-to- 
cold  relaxation  process  as  vibrationally  cold  N2  is  no  longer  being  provided  on  the  high  speed 
side.  This  effect  is  illustrated  above  in  Fig.  15b.  While  the  rotational  temperatures  for  the  two 
different  jet  compositions  match  well,  the  measured  non-equilibrium  is  much  stronger  when 
there  is  no  longer  N2  present  in  the  jet.  Measurements  with  the  H2/Ar  jet  are  particularly  difficult 
due  to  the  lack  of  N2  from  the  high-speed  side  of  the  shear  layer  to  provide  scattered  signal;  this 
is  why  the  profile  for  these  runs  does  not  extend  as  far  into  the  jet. 

Values  for  freestream  density  and  mass  entrainment  ratios  are  presented  below  in  Table  4.  If 
more  mass  from  the  cold  jet  is  entrained  (Em  >  1),  the  mixing  layer  temperature  profile  will  be 
depressed  such  that  the  hot  side  will  be  further  from  the  mean  temperature.  The  hot-to-cold 
process  dominates  entrainment  in  all  cases,  which  agrees  with  the  behavior  seen  in  the 
experimental  shear  layer  profiles  as  well  as  the  CFD  results. 14  However,  the  cold-to-hot  process 
is  still  present  and  should  influence  the  mean  vibrational  temperature  profiles  by  reducing  the 
apparent  non-equilibrium  measured,  as  seen  in  Fig.  15b. 

Table  4:  Calculated  shear  layer  density  and  mass  entrainment  ratios 


Case 

S 

h2/n2 

1.43 

1.67 

H2/Ar 

1.31 

1.75 

H2/CH4 

2.14 

1.37 

The  temperature  profile  at  a  location  farther  downstream  is  presented  in  Fig.  16.  This  location 
corresponds  to  the  farthest  downstream  location  where  the  flame  is  never  present.  At  this 
location,  the  jet  potential  core  is  just  collapsing  as  there  is  no  flat  section  in  the  temperature 
profile  corresponding  to  a  region  of  pure  jet  fluid;  therefore,  the  entire  jet  width  can  be  probed 
with  Raman  scattering  since  the  N2-rich  outer  flow  has  penetrated  to  the  center  of  the  jet.  The 


16 

Distribution  A  -  Approved  for  Public  Release 


shear  layer  profiles  on  each  side  of  the  jet  show  similar  behavior  to  the  upstream  profiles 
presented  above,  which  is  expected  since  this  location  is  just  beyond  the  shear  layer  region  of  the 
jet.  It  is  important  to  note  that  non-equilibrium  will  not  immediately  disappear  downstream  of 
this  point.  The  driving  force  behind  non-equilibrium  generation  is  the  temperature  gradient  in  the 
flow,  which  is  still  present  beyond  the  potential  core.  The  non-equilibrium  will  simply  decay  in 
magnitude  as  the  jet  centerline  temperature  rises  to  the  outer  flow  temperature  due  to 
downstream  turbulent  mixing.  This  behavior  is  also  observed  in  our  simulation  results  presented 
in  Ref.  14  (see  also  below)  where  it  is  shown  there  that  non-equilibrium  at  the  centerline  begins 
near  the  collapse  of  the  potential  core  and  continues  until  far  downstream  where  combustion 
products  have  been  entrained  to  the  jet  centerline. 

Temperatures  were  not  extracted  from  Raman  scattering  spectra  in  the  region  where  the  flame  is 
intermittently  present  due  to  rapid  turbulent  fluctuations.  Owing  to  the  intermittent  presence  of 
the  flame,  these  regions  have  a  widely  spread  bimodal  temperature  PDF,  where  the  two  peaks 
will  be  around  the  freestream  (1000  K)  and  post-flame  (2400  K).  As  seen  previously,  the  Raman 
fitting  procedure  will  only  converge  to  a  meaningful  temperature  value  if  the  temperature  PDF 
within  the  measurement  volume  is  sufficiently  narrow.  However,  measurements  can  be  taken 
downstream  of  this  region.  Results  at  such  a  downstream  station  are  plotted  in  Fig.  17.  The 
acquisition  of  spectra  at  this  axial  location  required  a  15  s  integration  time.  The  figure  at  left 
shows  a  mean  profile  through  both  sides  of  the  flame,  which  is  seen  to  be  symmetric,  as 
expected.  It  is  also  clear  that  the  non-equilibrium,  which  is  seen  in  the  non-reacting  regions 
upstream,  is  completely  gone  at  the  location  of  peak  temperature  in  the  radial  profile.  This  lack 
of  non-equilibrium  is  expected  since  a  high  concentration  of  flame  products  such  as  H20  should 
quickly  relax  the  non-equilibrium.  Non-equilibrium  is  detected  at  the  jet  centerline  just  as  it  is  in 
the  upstream  profile  of  Fig.  16.  Our  simulation  results  show  a  similar  non-equilibrium  zone  at 
the  centerline  near  the  height  of  the  mean  flame  base.  Figure  17b  shows  a  zoomed- in  look  at  the 
peak  temperature  location.  The  peak  measured  temperature  of  1835  K  is  240  K  less  than  the 
adiabatic  flame  temperature  predicted  for  this  fuel,  which  is  reasonable  considering  the 
measurement  is  time-averaged  and  the  flame  base  fluctuates  substantially. 


Figure  16:  Radial  mean  profiles  of  N2  rotational  and  vibrational  temperature  approximately  four  diameters 

downstream  of  the  nozzle  exit. 
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Figure  17:  Radial  trace  of  N2  temperatures  7  diameters  downstream  of  nozzle  exit  (downstream  of  flame  base)  (a) 
Profile  across  entire  flame,  (b)  zoomed-in  view  to  show  data  near  the  peak  temperature 


2.  Theoretical  effort 

Bimolecular  chemical  reactions  may  be  dissociation  or  exchange  type  and  are  usually  written  as 
AB  +  CD  — »  A  +  AB  +  CD  and  AB  +  CD  — >  AC  +  BD  respectively.  This  description  is  adequate 
to  describe  reaction  rates  under  thermal  equilibrium  conditions  (i.e.  Maxwellian  velocity 
distributions  and  Boltzmann  energy  distributions  for  internal  states).  It  must  be  extended  for 
conditions  of  thermal  non-equilibrium  and  for  dissociation  reactions,  for  example,  we  must  write 
AB(vj,  Jj)  +  CD(v2,./2)  — >  A  +  B  +  CD(v3,  J3). 


Because  the  reaction  cross-sections  are  state  dependent  the  overall  reaction  rate  must  be  obtained 
by  summing  over  all  product  internal  states  and  averaging  over  reactant  internal  states: 


Rf=  kfnABnCD  =  nABnCD  XXX  ~  ~ ~ ~ ~  ~  ~ ~ ~ ~  2  3>  d)Crel  ^ 


V3./3  F|  J\  VjJ 2 


n 


AB 


*CD 


Here  the  state-specific  reaction  cross-section  (dissociation  cross-section  in  this  example)  is 
averaged  over  the  relative  speed  distribution  function,  which  may  be  assumed  to  be  Maxwellian 
except  under  extreme  non-equilibrium  conditions: 

oo 

cr(l,  2  — »  3,  d)crel  ■  =  J  <7(1, 2  ->  3,  d)crel%M  (g;  T)dcrel . 

0 


Reaction  rates  in  non-equilibrium  flows  differ  from  equilibrium  rates  because  the  internal  states 
of  the  reactants  are  not  described  by  Boltzmann  distributions.  One  requires  state-specific 
reactions  rates  for  the  most  detailed  description,  but  this  is  computationally  intractable  for 
reacting  flows  in  complex  geometries.  Multi-temperature  models  of  reactions  are  an  intermediate 
description  which  assumes  that  internal  energy  modes,  in  particular  the  vibrational  modes,  are 
Boltzmann  distributed  with  vibrational  temperature  Tv,  and  Tv  *  Tr  =  T,.  Quasi-classical 
trajectory  (QCT)  calculations  on  potential  energy  surfaces  (PESs)  determined  from  ab  initio 
quantum  chemistry  must  be  used  to  determine  state-specific  rates.  Depending  on  the  level  of 
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detail  required  the  initial  states  of  the  reactants  are  sampled  over  appropriate  Boltzmann 
distributions  to  obtain  suitably  averaged  rates. 

QCT  code  development 

We  have  developed  and  tested  a  highly  scalable  parallel  algorithm  to  run  QCT  trajectories  on 
general  PESs.  The  program  is  designed  to  distribute  the  computational  effort  over  available  cores 
and  adjusts  the  sampling  automatically  to  match  a  desired  global  uncertainty  as  shown  in  Figure 
18.  Each  process  samples  many  trajectories  and  inter-process  communications  are  designed  as 
checkpoints.  It  should  be  noted  that  the  number  of  checkpoints  is  independent  of  the  number  of 
active  cores. 


Master  Node  Path  (  Worker  Node  Path 

Master  &  Worker  Node  Path  |  MPI  Communication  Path 


Figure  18.  Schematic  diagram  of  program  workflow 

The  program  was  tested  by  running  5  million  trajectories  for  the  N2+N2  system  and  demonstrated 
linear  scaling  of  the  computational  time  with  the  number  of  active  cores.15  Figure  19  shows  the 
results  of  two  test  runs  to  calculate  inelastic  cross-sections  for  N2-N2  collisions.  For  Test  A  there 
are  200  trajectories  simulated  in  each  iteration  i.e.  each  batch  of  trajectories.  In  Test  B,  there 
were  800  trajectories  per  iteration,  resulting  in  a  reduction  by  a  factor  of  four  in  the  number  of 
MPI  calls  because  there  are  fewer  instructions  from  the  master  node.  The  total  computational 
time  as  a  function  of  the  number  of  processors  is  shown  in  Fig.  19a.  Both  tests  show  excellent 
weak-scaling  as  the  number  of  processors  was  increased  from  24  to  4096,  with  -10%  increase 
over  the  minimum  for  Test  A,  and  -4%  increase  over  the  minimum  for  Test  B.  Figure  19b  shows 
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the  variation  of  the  average  CPU  time  per  trajectory  with  the  number  of  trajectories  simulated 
per  processor.  This  plot  shows  that  faster  processors  generally  simulate  more  trajectories,  so 
slow  processors  do  not  act  as  a  bottleneck  to  the  computations. 


ivproc  iraiecioi 

(a)  (b) 


Figure  19.  Demonstration  of  parallel  scaling  up  to  4096  cores  for  QCT  calculations  of  N2+N2  collisions 


Code  verification 

The  code  was  verified  by  simulating  500  million  trajectories  on  the  N4  PES  of  Paukku  et  al.16 
and  compared  to  the  results  of  Bender  et  al.17  who  used  the  same  PES  and  assumed  an  internal 
(ro-vibrational)  temperature  Trv  that  is  possibly  different  from  the  translational  temperature  Tt.  In 
our  calculations  results  were  obtained  on  a  relatively  coarse  grid  Tt  and  Trv  and  interpolated  onto 
a  full  surface  using  the  interpolation  procedure  of  Wang,  et  al.18  The  QCT  calculations  were  run 
with  Monte  Carlo  sampling  over  initial  conditions  (orientation  of  molecules,  vibrational  phase, 
etc.).  Details  are  given  in  the  doctoral  dissertation19  of  Stephen  Voelkel,  one  of  the  graduate 
students  supported  by  this  grant.  For  an  initial  state  E,  and  desired  final  state  E,\  with  Nr 
trajectories  resulting  in  the  desired  transition  out  of  N  sampled  trajectories,  then  the  probability 
of  the  transition  P{E,  — »  £')  »Nr/N,  and  the  precision  of  the  estimate  improves  for  large  sample 
numbers.  The  uncertainty  can  be  expressed  quantitatively  as: 


sp(t->e)= 2 


nJ  n-n, 
jv(jv„(jv-i)J 


Selective  sampling  for  increased  efficiency 

We  have  developed  a  new  sampling  procedure  to  reduce  computational  cost  when  sampling  for 
processes  that  have  low  probability  of  occurring.  This  is  critical  to  computing  reaction  rates  at 
modest  temperatures  (2,000-6,000  K)  which  are  commonly  found  in  engineering  applications. 
For  low  probability  outcomes,  the  standard  Monte-Carlo  technique  requires  a  very  large  number 
of  samples  in  order  to  obtain  reasonable  statistics  and  corresponding  precision  on  desired 
outcome.  Low  probability  reactive  events  correspond  to  high  activation  energy  processes  which 
require  relative  translational  energies  that  are  only  in  the  high  energy  tail  of  the  Maxwellian 
distribution.  We  have  implemented  a  scheme  that  only  samples  trajectories  beyond  the  estimated 
cut-off,  and  accounts  analytically  for  all  the  trajectories  that  are  below  the  required  threshold. 
Thus  the  effective  number  of  trajectories  can  be  several  orders  of  magnitude  larger  than  the 
number  of  simulated  trajectories  which  greatly  increases  efficiency. 
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Figures  20a, b  compare  the  measured  convergence  rate  of  the  relative  uncertainty  in  the  N2 
dissociation  rate  coefficient  for  standard  sampling  and  selective  sampling.  At  high  temperature 
(30,000K)  the  rate  coefficient  is  high  and  there  is  no  change  in  the  required  number  of 
trajectories  when  using  selective  sampling.  However,  at  6,000  K,  the  dissociation  rate  is  low,  and 
the  number  of  trajectories  required  to  get  ~10%  error  in  the  dissociation  rate  is  reduced  by  a 
factor  of  100  from  10  to  10  .  The  improvement  will  be  even  larger  at  lower  temperatures.  We 
believe  this  efficiency  improvement  will  be  very  useful  when  using  QCT  for  many  engineering 
applications. 


Figure  20.  Variation  in  estimated  error  of  computed  N2  dissociation  rate  with  number  of  trajectories  for  (a)  Standard 
sampling  and  (b)  Selective  sampling  at  T=  6,000,  13,000,  and  30,000  K. 


T  /  1000  K 
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Figure  21.  Verification  of  the  QCT  code  by  comparison  of  our  equilibrium  results  with  other  calculations  and 

experiment. 

Figure  21  compares  the  thermal  equilibrium  dissociation  rate  computed  via  QCT  using  selective 
sampling  in  the  range  6,000  K  <  T  <  60,000  K.  The  uncertainties  are  small  even  at  6,000  K 
because  of  the  selective  sampling  procedure.  Our  results  are  virtually  indistinguishable  from  the 
QCT  results  of  Bender,  et  al.17  in  the  range  8,000  K  <  T  <  30,000  K  obtained  using  the  same 
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potential  energy  surface.  Bender  does  not  use  selective  sampling  and  did  not  present  calculations 
at  temperatures  below  8,000  K. 


10  20  30  40  50  60 

t;/ioook 

(a) 


10  20  30  40  50  60 

rt/1000K 

(b) 


10  20  30  40  50  60 

r*/1000K 
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Figure  22.  (a)  Non-equilibrium  N2  dissociation  rate  kd(Tt  =  Tr,  Tv)  computed  from  QCT  code  (b)  Non-equilibrium 
rate  from  Park  model,  (c)  kd{Tt,  Tr  =  Tv )  computed  using  the  QCT  code. 


Figures  22  (a)  and  (b)  compare  the  N2  non-equilibrium  dissociation  rate  in  the  range  6,000  K  < 

90 

Th  Tv  <  60,000  K  computed  via  QCT  and  compares  the  results  with  Park’s  model 
Kark=Kq  ( Teff )  where  Te(f.  =^TtTv  .  We  use  Tr  =  T,  for  the  results  shown.  Note  the  logarithmic 

scale  on  the  plots.  The  two  models  agree  along  the  diagonal  Tv=Tt  by  design.  It  can  be  seen  that 

Park’s  model  works  reasonably  well  at  high  temperatures  but  there  are  significant  discrepancies 
at  low  temperatures.  Figure  22  (c)  shows  the  non-equilibrium  dissociation  rate  assuming  Tr  =  Tv. 
The  QCT  code  has  been  used  to  compute  the  non-equilibrium  dissociation  rate  kd(Tt,  Tr,  Tv)  for 
N2  and  an  approximate  three  temperature  rate  model  is  proposed  that  can  be  used  in  CFD  codes 
that  currently  use  Park’s  model. 


K  (T„Tr,T,)  =  K"  s  crc;  c,  +c2+c,=l. 


The  accuracy  of  the  model  is  quantified  using  the  root  mean  square  (RMS)  of  the  relative  error 
of  the  fit,  i.e.  e  =  ^k^q  {Te)-kd{Tt,Tr,Tv)^lkd{Tt,Tr,Tv^  .  The  fit  coefficients  for  several 

temperature  ranges  are  given  in  Table  5.  A  single  fit  over  the  entire  temperature  range  gives 
results  to  no  better  than  a  factor  of  2  (nearly  100%  error),  but  much  higher  accuracy  is  possible 
for  restricted  temperature  ranges  that  are  relevant  to  re-entry  flows. 


Table  5:  Coefficients  for  effective  temperature  model  over  selected  temperature  ranges 


T range  [K] 

Cl 

C2 

C3 

RMS  error  (%) 

6,000-60,000 

0.3133 

0.2787 

0.4080 

98.6 

13,000-40,000 

0.3092 

0.2258 

0.4650 

22.0 

30,000-60,000 

0.4692 

0.1777 

0.3531 

2.86 

Hydrogen  oxidation  reactions 

The  QCT  code  was  applied  to  calculate  non-equilibrium  reaction  rates  for  three  key  reactions  in 
the  hydrogen  combustion  mechanism: 
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H  +  02  — »  O  +  OH;  O  +  H2  — »  H  +  OH;  H2  +  OH  — »  H  +  H20 

2 1 

The  potential  energy  surfaces  and  other  computational  details  are  presented  in  Voelkel,  et  al. 
Non-equilibrium  reaction  rate  coefficients  kr{  T,  Tv)  for  these  reactions  were  computed  assuming 
Tr=  T,  =  T.  For  reaction  3  (H2  +  OH)  we  examined  the  effect  of  vibrational  non-equilibrium  of 
H2  and  OH  independently.  In  order  to  apply  the  results  conveniently  in  CFD  codes  and  efficiency 
function  is  defined  (  7\  T  )  =  kr  ( T,  Tv ) j kerq  (T).  An  important  conclusion  of  this  study  is  that 

exchange  reactions  of  the  type  studied,  are  much  less  sensitive  to  vibrational  non-equilibrium 
than  is  predicted  by  the  standard  Park  model  that  has  often  been  used  in  absence  of  better  data. 
Further  the  efficiency  functions  are  reaction  and  species  dependent.  This  can  be  seen  in  in  Fig. 
23  where  the  efficiency  function  contours  are  plotted  as  a  function  of  T,,/T  and  T  for  the  reactions 
above. 
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Figure  23.  Contour  maps  of  the  efficiency  function  showing  dependence  on  Tv/T  and  T  for  key  hydrogen  combustion 

reactions. 


While  there  is  generally  an  enhancement  in  the  reaction  rate  ( <p  >  l)  when  Tv  >  T9  the  magnitude 

depends  on  the  reaction.  The  first  two  reactions  are  more  sensitive  to  vibrational  non-equilibrium 
at  low  T,  while  the  H2  +  OH  reaction  is  more  sensitive  to  vibrational  non-equilibrium  of  H2  at 
high  T,  and  is  very  insensitive  to  vibrational  non-equilibrium  of  OH. 

Figure  24  shows  the  efficiency  surfaces  for  these  reactions  as  a  function  of  T  and  Tv  in  red.  The 
QCT  results  are  compared  to  the  efficiency  that  would  be  obtained  using  the  standard  Park 
model  (blue  surfaces)  and  the  coupled  vibration  chemistry  vibration  (CVCV)  model  (aqua 
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surfaces)  that  was  developed  by  Knab  et  al.  ’  to  describe  non-equilibrium  chemistry  for 
hypersonic  flows.  In  order  to  facilitate  use  of  these  results  in  standard  CFD  codes,  we  fit  the 

QCT  surface  to  an  “optimized”  Park-type  model,  <p(T,Tv  j  =  keq  ( Teff  )  j keq  (  71 ) ;  Teff 

The  standard  Park  model  corresponds  to  C=  1.  Neither  the  CYCY  model  nor  the  standard  Park 
model  matches  the  QCT  calculations  well.  The  optimum  value  of  C  for  these  exchange  reactions 
is  less  than  0.2,  and  the  particularly  low  value  Copt  for  the  Th  +  OH(rv)  reaction  reflects  how 
insensitive  the  rate  coefficient  is  to  vibrational  non-equilibrium  of  OH.  Table  6  shows  the  RMS 
errors  of  the  model  efficiencies  surfaces  relative  to  the  QCT  values. 


H  +  02(rv)^  O  +  OH;  ^  =  0.182 


0  +  H2(Tv)— >•  H  +  OH;  Copt  =  0.080 


H2  (Tv )  +  OH  H  +  H20 ;  CPt  =  0.097 


H2+OH(rv)^  H  +  H20 ;  Copt  =  0.013 


Figure  24.  Efficiency  function  surfaces  as  a  function  Tv/T  and  T  for  key  hydrogen  combustion  reactions.  Red  -  QCT 
computation;  Cyan  -  CVCV  model;  Blue  -  Standard  Park  model;  Yellow  -  Optimized  Park-type  fit  to  QCT 

calculated  surface. 
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Table  6:  RMS  error  of  the  modeled  efficiency  functions  compared  to  the  QCT  results 


T  range  [K] 

Copt 

This  work 
k(T,Tv)  =  keq(Teff) 

Teff  =  {tT^ 

Original  Park 

CVCV 

H  +  O  2{Tv) 

0.182 

0.012 

0.211 

0.283 

0  +  H2(7V) 

0.080 

0.010 

0.435 

0.108 

OH(7)  +  H  2(Tv) 

0.097 

0.016 

0.242 

0.047 

OH(Tv)  +  H2(7) 

0.013 

0.003 

0.306 

0.107 

In  addition  to  the  purely  empirical  two -temperature  fit  described  above,  we  have  generated  a 
more  physics  based  species-specific  efficiency  function  similar  to  the  model  presented  by 
Aresntiev  et  al.24  This  is  based  on  the  forward  and  reverse  activation  energy  of  the  reaction  as 
well  as  the  initial  and  final  vibrational  energy  of  the  products  and  reactants,  respectively.  For  a 
specific  reaction,  we  define  the  analytical  state-dependent  efficiency  function: 


^(v,v';T)  =  c1(r)ex p 


K  ( v )  -  max  ( s~  ( v' )  -  s~ ,  0 ) 


kT 


where  £v+  (v)  and  sv  (V)  are  the  vibrational  energy  of  the  reactants  and  products,  respectively, 
s~a  is  the  activation  energy  of  the  reverse  reaction,  k  is  the  Boltzmann  constant,  c\(T)  is  chosen 
so  that  k(T;  Tv  =  T)  =  ke(T),  and  c2(7)  is  a  free  parameter  set  to  (OA-T/IO4  j  based  on  QCT 
calculations. 


From  the  state-specific  function  (p(v,v';T'j  we  can  calculate  suitably  summed  and  averaged 
efficiency  functions,  cp ( v; T)  =  ^ cp ( v, v'; T) and cp (T, Tv )  =  ^/v (v, Tv ) cp ( v; Tv ), where  fv ( v,  Tv ) is 

v'  v 

the  fractional  population  in  state  v  computed  assuming  a  Boltzmann  distribution  at  Tv. 

The  state-specific  model  (p(v; T )  was  compared  to  (p{v\T )  determined  directly  from  QCT 

calculations  for  the  reaction  H  +  02(v  =  0;  1;  2;  3;  4;  5)  — >  O  +  OH  using  the  PES  developed  by 
Melius  and  Blint25  as  distributed  on  POTLIB.26  Approximately  14  million  trajectories  were 
simulated,  with  more  trajectories  simulated  at  low-probability  states  to  minimize  the  uncertainty 
(defined  as  two  standard  deviations  from  the  mean).  Figure  25  shows  the  proposed  model 
compared  to  the  QCT  results.  The  proposed  model  closely  matches  the  QCT  results,  especially 
as  T  increases.  The  model  is  not  a  good  match  for  high  v  at  low  T.  However,  the  population 
fraction  fv  of  high-lying  v  is  very  small  at  low  T  so  the  departure  from  the  model  fit  has  little 
effect  on  the  temperature  dependent  efficiency  function. 

The  temperature-dependent  model  <p(T,Tv)  was  compared  to  the  standard  T-Tv  model  of  Park20 

and  the  CVCV  model23  for  two  of  the  reactions  in  the  hydrogen-air  combustion  mechanism: 
H+02(rv)  — »  O+OH  and  O  +  H2(TV)  — »  H  +  OH.  Additionally,  these  models  were  compared  to 
<p(T,Tv) directly  calculated  via  QCT  as  described  above.  The  resulting  plots  are  very  similar  to 
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those  shown  in  Fig.  24  above.  For  both  reactions,  the  QCT  results  show  that  vibrational 
nonequilibrium  has  a  smaller  effect  than  the  Park’s  model  and  the  CYCY  model  for  both  Tv>  T 
and  Tv  <  T.  Overall,  the  physics-based  model  <p(T,T  "j  matches  the  QCT  results  better  than  the 

Park  model  and  the  CVCV  model  except  when  Tv  >  T  and  T  <  1400  K  for  the  reaction  H  + 
O2 (Tv)  — »  O  +  OH,  where  Park’s  model  is  closer  to  the  QCT  results. 


Figure  25.  Vibrational  energy  exchange  probability  as  a  function  of  T  and  Tv()2  or  Tv^2. 

We  also  made  QCT  calculations  of  the  probability  of  V-V  transfer  rate  between  N2  and  O2  using 
the  N2O2  potential  energy  surface  from  POTLIB.  These  calculations  were  motivated  by  the  jet 
flame  experiments  described  in  Section  1,  which  showed  that  in  the  turbulent  mixing  layer 
upstream  of  the  flame  base,  the  vibrational  temperature  of  oxygen  equilibrated  rapidly  with  the 
translation-rotational  temperature,  but  the  vibrational  temperature  of  nitrogen  remained  out  of 
equilibrium  (see  Fig.  15a  above).  Simulations  of  this  flow  using  Pw  =  10  2  (the  value  of 
vibrational  energy  exchange  probability  from  the  standard  CVCV  model)  did  not  reproduce  the 
experimental  data.  Simulations  (described  in  Section  3  below)  were  run  with  smaller  values  of 
Pw  and  best  agreement  was  found  for  Pw=  10  4.  The  results  of  the  QCT  calculations  shown  in 
Fig.  26  confirm  that  in  the  temperature  range  of  the  experiment  (500  -  1000  K)  the  vibrational 
exchange  probabilities  are  indeed  of  the  order  10  4. 


Figure  26.  Vibrational  energy  exchange  probability  as  a  function  of  T  and  Tv  0l  or  Tv  Kl. 


26 

Distribution  A  -  Approved  for  Public  Release 


3.  Computational  effort 

Computations  for  this  project  involved  both  LES  and  DNS  flow  simulations  of  the  experimental 
configuration,  a  model  of  a  scramjet  combustor  (Hyshot-II),  and  also  of  detonation  waves.  We 
have  reported  on  a  number  of  simulations  in  interim  reports  which  have  appeared  in  several 
publications.  Here  we  report  on  simulations  performed  in  the  final  year  of  the  grant. 

Simulations  of  the  experimental  set-up 

Detailed  computations  of  a  fuel  jet  were  undertaken  to  support  the  experimental  effort.  The 
computational  domain  is  shown  schematically  in  Figure  27.  The  fuel  nozzle  diameter  ( D )  was  8 
mm,  the  coflow  diameter  was  100  mm  to  match  the  experiment.  The  flow  Reynolds  number  Red 
=  32,000  and  the  axisymmetric  LES  spatial  grid  used  was  (Nx,Nr,Ng)  =  (392,320,16). 

Fuel  jet 
Ma  =  0.568 
p  =  0.221  kg/m3 
T0  =  500  K 
Vjet  =  1000  m/s 
Th2  =  0.1765 
7n2  =  0.8235 

Air  coflow 
Ma  =  0.05 
p=  0.353  kg/m3 
Ta=  1000  K 

E yjfiow  30  m/s 

7o2  =  0.233 
7n2  =  0.767 


Figure  27.  Schematic  diagram  of  the  geometry  and  flow  parameters  used  in  LES  flow  computations 

The  UTComp  compressible  flow  solver  uses  a  fifth-order  WENO  scheme  with  characteristic 
reconstruction  to  compute  the  nonlinear  momentum  terms,  while  a  fourth-order  central  scheme  is 
used  for  the  diffusion  terms. 27,28’29’30’3 1,32  The  computational  time  step  was  50  ns.  The  vibrational 
temperature  of  each  species  is  tracked  independently.  The  standard  Millikan  and  White 
correlation33  is  used  to  model  V-T  relaxation.  The  chemistry  was  modelled  using  the  9-species, 
19-step  mechanism  of  Mueller,  et  al.34  We  use  the  QCT  derived  rates  for  the  three  key  reactions 
as  described  in  the  previous  section.  For  the  remaining  16  reactions  in  the  mechanism  we  use  the 
CVCV  model  of  Knab,  et  al.23  We  also  use  the  QCT  derived  rates  for  N2  dissociation,  though 
this  reaction  is  not  significant  for  these  flow  conditions.  The  simulation  was  run  for  8  hours  on 
2000  cores  with  a  CFL  number  of  0.9.  Figure  28  shows  sample  results  for  instantaneous 
translational  temperature,  and  N2  vibrational  temperature  obtained  from  the  LES  computations. 
The  computed  flow  is  more  lifted  than  observed  experimentally  because  the  co-flow  velocity 
used  in  the  simulation  (30  m/s)  was  higher  than  the  experimental  value  of  ~  1  m/s.  Stable 
numerical  solutions  could  not  be  obtained  for  a  very  low  co-flow  velocity. 
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The  computations  could  not  reproduce  the  different  profiles  of  mean  vibrational  temperatures  of 
N2  and  O2  observed  experimentally  (See  Fig.  15a  above)  when  using  Pw =  10“2  for  the  CVCV 
model  as  in  Knab,  et  al.  Because  the  Knab  model  was  developed  for  reentry  flows  at  much 
higher  temperatures,  we  varied  the  value  of  Pw  until  a  good  match  was  obtained  between  the 
computed  mean  profiles  and  the  experimental  data.  Because  of  the  difference  in  ffeestream 
velocity  between  experiment  and  simulation  we  compare  mean  profiles  by  normalizing 
downstream  distance  by  the  potential  core  length,  and  the  data  are  presented  in  scaled  variables: 
temperatures  are  normalized  by  the  temperature  difference  between  jet  and  coflow  ffeestreams, 
and  radial  distances  are  normalized  by  the  shear  layer  width.  Figure  29  shows  the  computed 
mean  profiles  for  several  values  of  Pw  and  superposed  on  the  experimental  data  measured  2  jet 
diameters  downstream  of  the  nozzle  exit.  As  expected,  the  difference  between  the  vibrational 
temperature  profiles  of  O2  and  N2  increases  as  Pw  decreases  until  Pw  is  less  than  10”4,  with  little 
change  thereafter.  Figure  30  shows  profiles  of  rVjN2  and  rr>N2  measured  4  diameters  downstream 
of  the  jet  exit,  just  upstream  of  the  flame  base  (see  Fig.  16).  The  computed  mean  profiles  of  T 
and  rv,N2  for  this  location  are  superposed  on  the  data  for  Pw  =  10  2  and  10  4.  Again,  the 
agreement  between  data  and  experiment  is  greatly  improved  when  Pw  is  reduced  to  10”4.  This 
comparison  between  computations  and  experiment  and  experiment  motivated  the  QCT  study  of 
W  transfer  between  N2and  O2  described  in  the  previous  section.  As  noted  above  (See  Fig.  26) 
the  QCT  calculations  indicate  that  Pw  <  10”4  for  these  flow  conditions. 


T  [K]  Tv,  N2[K] 
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Figure  28.  Sample  results  from  LES:  Instantaneous  (left)  translational  temperature  T;  and  (right)  N2  vibrational 

temperature,  Tv  K-, 


Figure  29.  Comparison  of  experimentally  measured  profiles  of  mean  Tr  Kr  Tro2,  Tvn2,  Tv,o2  2  jet  diameters 
downstream  of  the  nozzle  exit  with  computed  mean  values  at  the  same  location  for  different  values  of  Pvv. 
Rotational  temperatures  are  in  equilibrium  with  each  other  and  with  translational  temperature  measured  by  Rayleigh 
scattering.  Good  agreement  for  Pyy  =  KT4  with  little  change  in  the  profiles  if  the  value  is  reduced  further. 

Detailed  computations  were  performed  on  the  Hyshot  II  configuration  to  implement  the  non¬ 
equilibrium  reaction  rates  determined  from  the  QCT  calculations  in  a  more  realistic  scramjet 
configuration  mimicking  high  altitude  flight  at  Mach  8.  Figure  31  shows  the  geometry  consisting 
of  an  18°  ramp  leading  to  the  isolator  and  combustor.  Two  different  cases,  with  and  without 
thermal  nonequilibrium,  are  simulated.  The  computational  domain  is  divided  into  two  parts,  with 
the  initial  bow  shock  created  by  the  wedge  simulated  using  a  standard  RANS  solver.  For  the 
equilibrium  condition,  the  solver  directly  provides  the  post-shock  conditions  based  on  a  shock 
capturing  scheme.  For  the  nonequilibrium  case,  the  vibrational  temperature  through  this  shock  is 
frozen  while  the  translational  temperature  is  allowed  to  increase  in  order  to  preserve  total  energy. 
The  computational  domain  starts  further  downstream  at  x  =  0  cm,  and  extends  to  x  =  36.5  cm. 
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The  far-field  pressure  is  1197  Pa  with  air  density  of  0.014  kg/m  ,  while  the  far- field  temperature 
is  set  to  300  K.  The  flow  over  the  wedge  is  first  computed  using  a  RANS  model  to  estimate  the 
turbulent  boundary  layer  height  and  momentum  thickness  over  the  wedge  section.  These  values 
are  used  to  create  inflow  conditions  for  the  main  simulation  domain  based  on  a  separate  DNS  of 
a  periodic  turbulent  boundary  layer.  The  inflow  data  is  collected  from  this  auxiliary  DNS 
calculation  and  fed  to  the  actual  simulation  domain  used  in  this  work.  The  Reynolds  number 
based  on  the  inflow  boundary  layer  height  is  18,800  at  Mach  4.2  (post-shock).  For  the 
nonequilibrium  case,  the  flow  behind  the  first  bow  shock  is  modified  to  incorporate  the 
differences  in  vibrational  and  translational  temperature  using  the  ID  conservation  equations  with 
V-T  relaxation.  In  this  way,  the  inflow  to  the  computational  domain  contains  the  same  internal 
energy  for  both  cases.  A  bleed  section  is  introduced,  similar  to  that  in  the  experimental 
configuration,  to  naturally  remove  any  additional  shock  reflections  in  the  combustor  section.  The 
bleed  is  located  in  such  a  way  that  for  the  flow  Mach  number  considered,  the  foot  of  the  reflected 
bow  shock  is  captured  by  the  bleed  inlet.  The  RANS  simulation  for  the  wedge  section  is 
performed  using  a  realizable  k-s  model  on  FLUENT  running  on  8  processors  over  a  2D 
orthogonal  grid  of  ( nx ,  ny)  =  (2048,  512)  until  convergence  is  reached  (scaled  residuals  less  than 
0.001).  The  high-resolution  calculations  are  performed  using  the  compressible  flow  solver 
UTComp  described  briefly  above.  For  these  calculations  a  CFL  number  of  0.85  is  used  for  a  time 
step  of  3.9  ns. 

The  computational  domain  has  over  33  million  nodes  with  a  wall- normal  refined  grid  of 
(nx,ny,nz)  =  (2048,128,128).  The  walls  are  assumed  isothermal  at  300  K,  with  periodic  boundary 
conditions  in  the  spanwise  (z)  direction.  The  grid  size  is  Az  =  1 6  v ' .  In  the  ramp  section  Av  =  18 
y1  and  Ay  =  [l:17]y+,  with  y 1  estimated  to  be  ~10  pm.  A  preliminary  RANS  calculation  is 
performed  on  the  forward  portion  of  the  wedge,  and  the  turbulent  boundary  layer  on  the  wedge  is 
computed  from  an  auxiliary  DNS  computation  with  Re  =  18,800.  The  viscosity  is  determined 
using  Sutherland’s  law  while  thermal  diffusivity  is  obtained  from  a  constant  Prandtl  number  of 
0.7.  Both  cases  were  carried  out  on  4000  processors  using  MPI-based  parallelization  for  roughly 
120  wall-clock  hours.  The  centerline  averaged  flow-through  time  is  rc=  160  ps,  and  the  statistics 
were  sampled  over  4rc  after  reaching  statistical  stationarity.  We  confirmed  that  longer  averaging 
did  not  change  the  statistics  significantly. 

The  resolution  of  the  simulations  was  based  on  the  reaction  layer  thickness.  Note  that  for  the 
turbulent  flow  conditions,  a  DNS  will  require  another  order  magnitude  increase  in  computational 
mesh,  purely  based  on  isotropic  turbulence  scaling.  However,  this  geometry  has  two  differences. 
First,  the  reactions  are  sufficiently  slow  in  the  ignition  region  so  that  reaction  layers  are 
reasonably  broad.  For  instance  the  OH  layers  are  resolved  using  2-4  grid  points  for  the 
computational  grid  used,  even  at  highly  strained  conditions.  Second,  the  main  turbulence 
generation  is  through  the  interaction  of  the  lower  velocity  fuel  jet  and  the  crossflow.  Here,  the 
evolution  of  the  mixing/shear  layer  determines  the  small-scale  structure.  In  such  jet-in-crossflow 
problems,  the  fully-developed  turbulence  structure  is  not  established  until  further  downstream. 
Hence,  the  small  scales  are  bound  to  be  larger  compared  to  isotropic  turbulence  scaling.  This  is 
further  verified  by  the  size  of  the  reaction  zones.  Further  details  are  given  in  Fievet,  et  al.35,36 

The  simulations  show  that  vibrational  non-equilibrium  persists  to  the  combustor  inlet.  Excited 
vibrational  states  are  underpopulated  and  in  consequence  the  translational  temperature  is  higher 
than  in  an  equilibrium  flow.  Figure  32  shows  the  mean  temperature  profiles  on  the  centerline, 
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showing  that  the  translational  temperature  in  the  non-equilibrium  flow  is  300  K  higher  than  if  the 
flow  were  in  equilibrium  (1600  K  instead  of  1300  K). 


Figure  30.  Comparison  of  experimentally  measured  profiles  of  mean  T,  Tv,N2  4  jet  diameters  downstream  of  the 
nozzle  exit  with  computed  mean  values  at  the  same  location  for  different  values  of  Pw  The  measured  Tr  ^2  matches 
well  with  computed  translational  temperature(and  temperature  measured  by  Rayleigh  scattering).  Computed  mean 
profile  of  Ty^2  does  not  agree  with  experiment  for  Pw  =  10-4;  better  agreement  for  Pw  =  10“4. 
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Poo=l  197Pa 
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combustor  : 
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Figure  31.  Hyshot  II  geometry  being  simulated 


x=36.5 


The  effect  of  this  elevated  translational  temperature  on  the  three  hydrogen  reaction  studied  above 
is  shown  in  Table  7,  and  it  can  be  seen  that  the  equilibrium  rate  coefficients  for  reactions  2  and  3 
are  significantly  enhanced  at  the  higher  temperature.  The  nonequilibrium  rate  is  given  by 

k[T,Tvi)  =  (p(T,Tvi)keq  (r). Figure  32  shows  that  (p  is  not  very  sensitive  (<15%)  to  variations  in 

Tv  for  exchange  reactions  at  these  conditions.  Hence  it  does  not  completely  offset  the 
augmentation  of  the  reaction  rate  coefficient  by  the  higher  translational  temperature.  Thus,  for 
these  flow  conditions,  we  have  the  surprising  and  non-intuitive  result  that  the  key  reaction  rates 
are  higher  in  the  non-equilibrium  flow  than  in  an  equilibrium  flow.  Figure  33  shows  that  the  heat 
release  occurs  further  upstream  in  the  non-equilibrium  flow  and  the  flame  is  stabilized  further 
upstream. 
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Table  7:  Comparison  of  the  equilibrium  rates  at  1600  and  1300  K 


# 

Reaction 

keq  (1600  K)/*„  (BOOK) 

1 

H  +  02  ->  O  +  OH 

~1 

2 

O  +  H2  ->  H  +  OH 

1.23 

3 

OH+  H2  ->  H  +  H20 

1.45 

TJT 


Figure  32.  Variation  of  the  efficiency  factor  for  reactions  1-3  (Table  7)  at  1200  K. 


1G+10 
8E+09 
6E+09 
4E+09 
2B+09 
0 

0.1  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.2  0.21  0.22  0.23  0.24  0.25 

x[m) 

Figure  33.  Instantaneous  heat  release  rate  (W/m3)  computed  for  (top)  equilibrium  flow,  and  (bottom)  non¬ 
equilibrium  flow  simulations  of  the  Hyshot  II  combustor.  Hydrogen  injection  is  at  x  =  0.123  m  (see  Fig.  31). 

Since  each  species  is  described  by  a  specific  vibrational  temperature,  it  is  interesting  to  see  the 
evolution  of  these  vibrational  temperatures  as  the  flame  stabilizes.  Figure  34  shows  that  there  are 
significant  differences  between  the  relaxation  rates  of  the  vibrational  temperatures  of  individual 
species.  O2  exhibits  nonequilibrium  at  the  jet  inlet  due  to  the  inflow  conditions,  but  is  quickly 
equilibrated  when  the  initial  chemical  reactions  that  produce  H2O  occur.  This  is  similar  for  H2, 
which  exhibits  significant  nonequilibrium  during  the  expansion  process  close  to  the  injection 
point.  The  vibrational  temperature  of  H2O  shows  only  a  small  region  of  high  vibrational 
temperatures,  before  quenching  to  near-equilibrium  values.  It  is  important  to  note  that  the  flow 
away  from  the  reaction  zone  still  contains  substantial  nonequilibrium  due  to  the  presence  of  N2. 
Hence,  collisions  with  N2  molecules  transfer  nonequilibrium  to  other  species.  Further  details  are 
given  in  Ref.  35. 
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7(H2)  7(H20)  7(02) 

Figure  34.  Isocontours  of  mass  fraction  Yx=  0.1  for  (left)  H2,  (middle)  H20,  and  (right)  02.  In  each  case  the  mass 
fraction  contour  is  colored  by  the  ratio  of  species  vibrational  temperature  to  gas  translational  temperature.  The  ratio 
TVfX/T  for  each  species  on  the  symmetry  plane  is  also  shown.  Note  that  the  low  ratios  for  H2  and  H20  in  the  inlet  and 
isolator  far  upstream  of  the  injection  point  are  not  meaningful  as  there  is  no  significant  concentration  of  these 

species  far  upstream. 


4.  Summary  and  Conclusions 

There  has  been  a  fruitful  interaction  between  experiment,  theory,  model  development,  and 
numerical  simulations  under  this  grant.  High  dispersion  time-averaged  Raman  scattering 
measurements  were  used  to  identify  persistent  vibrational  non-equilibrium  of  nitrogen  in  high¬ 
speed  shear  layers  with  mismatched  gas  compositions.  Mixing-induced  thermal  non-equilibrium 
was  detected  in  the  shear  layer  upstream  of  the  turbulent  hydrogen  flame  in  N2  but  not  O2. 
Rotational  temperatures  of  the  two  species  agree  to  within  the  measurement  precision.  As 
expected,  the  non-equilibrium  is  relaxed  immediately  beyond  the  average  flame-base  location 
due  to  the  presence  of  combustion  products.  The  presence  of  non-equilibrium  is  confirmed  using 
Rayleigh  thermometry  images  to  quantify  the  effect  of  translational  temperature  variation  in  the 
Raman  measurement  volume.  These  results  were  compared  to  large-eddy  simulations  with 
vibrational  relaxation  effects  to  study  the  effect  of  interspecies  vibrational  energy  transfer 
models.  The  measurements  and  average  simulated  temperature  fields  could  only  be  made  to 
agree  when  the  interspecies  vibrational  coupling  is  very  weak.  This  required  a  two  order  of 
magnitude  change  in  the  probability  of  VV  transfer  that  was  widely  used  in  simulations  of  non¬ 
equilibrium  flow. 

A  highly  efficient  QCT  code  was  developed  under  this  grant  and  was  used  to  determine  non¬ 
equilibrium  reaction  rate  coefficients  for  three  key  reactions  in  the  hydrogen  combustion 
mechanism.  It  was  also  used  to  verify  the  experimentally  observed  slow  V-V  transfer  between 
nitrogen  and  oxygen,  which  allows  for  persistent  vibrational  non-equilibrium  prior  to  production 
of  water  with  consequent  rapid  vibrational  relaxation.  Simple  empirical  fits  were  used  to 
implement  these  QCT  based  non-equilibrium  rates  into  numerical  flow  simulations.  An 
improved  physics-based  model  was  also  developed,  that  shows  promise,  but  requires  further 
development. 
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Simulations  of  the  flame  from  a  high-speed  jet  in  a  heated  co-flow  showed  that  vibrational  non¬ 
equilibrium  delayed  flame  ignition.  However  for  hydrogen  injected  in  cross-flow  in  a  shock 
heated  air  stream  as  in  a  scramjet,  the  persistent  nitrogen  non-equilibrium  leaves  the  air  at  an 
elevated  translational  temperature  and  was  found  to  move  ignition  upstream  resulting  in  a  more 
stable  flame.  This  demonstrates  that  one  can  devise  geometries  where  non-equilibrium  can  be 
used  to  advantage  to  stabilize  a  flame  in  a  scramjet. 

Key  insights  arising  from  this  work  relative  to  scramjet  design  are 

•  N2  vibrational  nonequilibrium  plays  a  significant  role  in  determining  ignition  delay  (and 
hence  flame  stability)  under  scramjet  combustor  conditions 

•  The  nonequilibrium  effect  can  either  advance  or  delay  ignition  relative  to  equilibrium  flow 

•  Effect  arises  principally  from  impact  of  vibrational  nonequilibrium  on  translational 
temperature,  with  consequent  effect  on  reaction  rates 

•  Direct  effect  of  nonequilibrium  Tv  on  rate  coefficients  is  small  at  scramjet  conditions  because 
the  key  reactions  are  exchange  reactions  rather  than  dissociation  reactions 

•  The  nonequilibrium  can  be  controlled  via  strength  of  shocks  and  expansions,  and  location  of 
fuel  injection  relative  to  these  flow  features 

5.  Future  Work 

The  experimental  work  in  the  high  speed  jet  should  be  extended  to  enable  single-pulse  collection 
of  Raman  vibrational  spectra.  A  multi-pass  laser  configuration,  with  a  higher  energy  laser  pulse 
should  enable  single  pulse  measurements  of  vibrational  temperature  using  high  dispersion 
Raman  spectroscopy.  It  is  unlikely  that  rotational  temperature  can  be  obtained  via  single  pulse 
measurements,  because  the  O-  and  S-  branches  that  show  resolved  rotational  structure  are  much 
weaker  than  the  Q-branch.  Additional  measurements  to  advance  the  fundamental  understanding 
of  non-equilibrium  flows  and  processes  will  require  a  facility  that  operates  continuously  so 
reliable  measurements  can  be  made  of  high-lying  states  with  very  small  populations.  An 
inductively  coupled  plasma  torch  is  an  ideal  facility  for  such  experiments  as  it  produces  a  steady 
high  temperature  flow  with  no  electrode  contamination. 

The  QCT  code  is  highly  parallel  and  can  produce  state-specific  and  suitable  state-average  cross- 
sections  and  rate  coefficients  from  an  available  potential  energy  surface  (PES).  Currently,  the 
bottleneck  is  production  of  the  PES  from  ab  initio  quantum  chemistry  calculations  and 
generation  of  suitable  fits  to  the  high-dimensional  PES.  The  solution  of  challenging  engineering 
problems  under  extreme  conditions  that  are  of  interest  to  the  Air  Force,  would  be  greatly 
facilitated  if  the  PES  could  be  adaptively  generated  during  the  course  of  the  QCT  calculations. 
Additionally,  it  would  be  very  useful  if  quantitative  metrics  could  be  developed  to  specify  how 
accurately  the  PES  must  be  determined,  and  in  what  regions  of  the  high  dimensional  space  it 
must  be  known  accurately.  Presently,  the  fits  to  the  PES  are  determined  by  experts,  and  the 
metrics  used  often  focus  on  reproducing  static,  equilibrium  properties,  and  there  is  no  simple 
way  of  relating  this  to  accuracy  in  chemical  dynamics  calculations.  A  research  effort  in  this  area 
is  needed  to  improve  the  state-of-the  art. 
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